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We present a general approach to greatly increase at little cost the efficiency 
of Monte Carlo algorithms. To each observable to be computed we associate a 
renormalized observable (improved estimator) having the same average but a differ- 
ent variance. By writing down the zero-variance condition a fundamental equation 
determining the optimal choice for the renormalized observable is derived (zero- 
variance principle for each observable separately). We show, with several examples 
including classical and quantum Monte Carlo calculations, that the method can be 
very powerful. 
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Since the pioneering work of Metropolis et al. [|l| Monte Carlo methods have been widely 
used in many areas of natural sciences. At the root of Monte Carlo methods lies a very 
efficient stochastic method for calculating many-dimensional integrals (or sums) written 
under the general form 

^ JadxTT(x)0(x) 

JsdxTT{x) 

where 0{x) is some arbitrary observable (real- valued function) defined on the configura- 
tion space 5* (continuous or discrete) and 7t{x) some probability distribution. In Monte 
Carlo methods the integrals are evaluated using a large but finite set of configurations 
{a;'^*^}j=i^7vdistributed according to tt and generated by a step-by-step stochastic procedure 
(Markov chain) 

1 ^ 

<0>=-Y.0[x^^^]+60, (2) 
1=1 

where SO is the statistical error associated with the finite statistics. For a large enough num- 
ber N of Monte Carlo steps, standard statistical arguments lead to the following expression 
of the error 

60 = k'-^ (3) 



where K is some positive constant proportional to the amount of correlation between con- 
figurations and cr(0) a measure of the fiuctuations of the observable 



o[0) = V< 02 > -< O >l (4) 

In this Letter it is shown that by introducing a suitably renormalized observable 0{x) 
the statistical error can be drastically reduced and even suppressed, thus defining a zero- 
variance principle for the Monte Carlo calculation of observables. To realize this, a trial 
operator H and a trial function tl}{x) are introduced (a trial matrix and a trial vector in the 
discrete case). The operator H is supposed to be Hermitian (in all practical applications, 
real symmetric) and is chosen such that 

dyH{x,y)^^) = Q. (5) 

On the other hand, the trial function il){x) is a rather arbitrary function which is simply 
supposed to be integrable. Now, the renormalized observable 0{x) associated with the 
observable 0{x) is defined as follows 

d{x) = 0{x) + IdyHj^m ^ (g) 

V7r(x) 

As a direct consequence of Eq. (|l|) and of the very definition of the Hermitian operator H, 
Eq. (^, we have the important property 
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< >=< o > 



(7) 



In other words, both quantities 0{x) and 0{x) can be used as estimators of the desired 
average. However, the statistical errors, which are controlled by cr(0) and cr(0), can be very 
different. The optimal choice for {H^ip) is obtained by imposing the renormalized function 
to be constant and equal to the exact average. This leads to the following fundamental 
equation 

j dyH{x,y)ij{y) = -[0{x)- < O >]^H^) ^ a{d) = (8) 

At this point it should be emphasized that the idea of using renormalized estimators 
for reducing the variance is not new. A number of applications have been performed using 
various "improved" estimators having a lower variance (See, e.g., 0], 0). The basic idea is 
to construct new estimators by integrating out some intermediate degrees of freedom and, 
therefore, removing the corresponding source of fluctuations. However, to the best of our 
knowledge, no general and systematic approach based on a zero-variance principle and valid 
for any type of Monte Carlo methods has been proposed so far. 

In this work the following strategy is proposed. First, a Hermitian operator H verifying 
(1) is chosen. Second, some approximate solution of Eq.(^ is searched for. The various 
parameters entering are then optimized by minimizing the fluctuations of the renormal- 
ized observable over a finite set of points distributed according to vr and obtained from a 
short Monte Carlo calculation. Finally, a standard much longer Monte Carlo simulation is 
performed using 0{x) instead of 0{x) as estimator. 

Choice of H. Clearly, a large variety of choices are possible for the trial operator H. For 
Monte Carlo algorithms satisfying the detailed balance condition (in practice, the vast ma- 
jority of MC schemes) a very natural choice is at our disposal. Denoting p{x y) the 
transition probability distribution defining the Monte Carlo dynamics, the detailed balance 
condition is written as tt{x)p{x —>?/)= T^{y)p{y —>■ x) for all pairs (x,y) in configuration 
space. A most natural operator to consider is 



H{x,y) 



\ ^ y) - 6{x - y)]. (9) 



From the detailed balance condition it follows that the operator H is symmetric, H{x, y) = 
H{y,x). The fundamental property (^ is verified since the sum-over-final-states for a tran- 
sition probability is equal to one. For continuous systems Schroedinger-type Hamiltonians 
can also be considered 

"--\±i^^y(^) (10) 

where V{x) is some local potential constructed to fulfill condition (^: 

1 JL d\fM 
2x n(x) i=i '-'•^i 
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where d is the number of degrees of freedom. Note that in Eq. ([T0| ) H is written using the 
standard quantum-mechanical notation for a local Hamiltonian in the x-space reahzation. 
Choice of if). Once the operator H has been chosen, the optimal choice for is the exact 
solution of the fundamental equation. Of course, in practice only approximate solutions are 
available. What particular form to choose for ip is very dependent on the problem at hand, 
on the type of observables considered, and also on the form chosen for the trial operator 
H . However, a most important point to be stressed is that the global normalization factor 
associated with ip is a, pertinent parameter of the trial function. Minimizing the fluctuations 
of the renormalized function a{0) with respect to it, we get 

^ 0(x) f dyH{x,y)i,(y) 



a(df = (T(Of 5 — . (12) 

^ / / dyH{x,y)i,{y) y ^ 



The correction to cr(0)^ being negative we obtain the important result that, whatever the 
choice made for the trial function (even the most unphysical one!), the optimization of the 
multiplicative factor always leads to a reduction of the statistical error. 

Our first application concerns the Monte Carlo calculation of the internal energy of the 
standard 2D-Ising model at various temperatures and linear sizes L = 5, 10, 20, and 25. 
The observable considered is the energy function given by E{S) = — Y.<i,j> SiSj (coupling 
constant J = 1, sum limited to nearest neighbors, and periodic boundary conditions). The 
probability distribution is vr(S') = exp[—PE{S)] with P = l/ksT. Here, S = {Si, ■ ■ ■ , Sn) 
with Si = ±1, and N = Lx L is the total number of spins. Simulations have been performed 
using a Swendsen-Wang type algorithm |^ (non-local updates of clusters of spins). To con- 
struct the trial operator H we have chosen to use the transition probability distribution of 
Monte Carlo algorithms with local updates ("Heat-Bath" -type algorithms). The probability 
of flipping the spin = ±1 at site i is given by 

^ ^ ^0SiSi _j_ g f3SiSi ^ ^ 

where e=l (no flip) or -1 (flip), and Si is the sum of neighboring spin values. With this 
choice and using Eq.(^) the fundamental equation (^) can be rewritten under the form 

N 

Y,p{S ^ TiS)[Q{S) - Q{T,S)] = E{S)- <E> 

i=l 



^{S) = QiS)^7TiS) (14) 

where the application Tj {i = 1,---,N) describes a flip at site i and is defined by 
Tj(S'i, ■ ■ ■ , Si, - ■ ■ , Sn) = {Si, ■ ■ ■ , —Si, ■ ■ ■ , Sn). At /? = (T = oo) the transition probability 
distribution becomes constant and the exact solution is easily found to be Q{S) = E{S)/2. 
For finite temperatures some approximate solution has to be found. Here we introduce for 
Q{S) a polynomial expansion up to the fourth-order in the variables X = Y^iLi Si (mag- 
netization) and Y = J2f=i g{SiSi) ("generalized energy", the usual energy being recovered 
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for g{x) = x). Precisely, we have chosen the form Q{S) = e Z]n+m<4 Cnm-^"^"* where 
Z = Y.iLih{SiSi). The set of variational parameters of ijj consists of all coefficients c„m 
of the polynomial plus the ten possible values of functions g and h. All coefficients have 
been optimized by minimizing the fluctuations of the renormalized energy a{E) defined by 
(I) and calculated from 2000 to 5000 different spin configurations S''-*-' drawn according to 
TT. Finally, the last step consists in performing a long Monte Carlo simulation to compute 
accurately the various quantities. The number of clusters built varies from 10^ (for the 
larger size) to 2.10^ (for the smaller size). Results are presented in Table |I[ Three different 
temperatures have been considered. T = 3 corresponds to the low-temperature regime, 
T = Tc = 4/ ln(v^ + 1) is the critical temperature for the infinite lattice, and T = 8 is in 
the high-temperature regime of the model. At T = 3, our representation is extremely good 
whatever the size of the lattice considered. The variance associated with the renormalized 
energy is drastically reduced with respect to the bare value and the gain in computational 
effort can be as great as ~ 360. Here, the gain in computational effort is defined as the ratio 
of the squared statistical eTTOTS,{SE/SE) . In other words, according to Eq.(|) it represents 
the factor by which it would be necessary to increase the number of Monte Carlo steps in 
the standard approach to get the same accuracy. Note that for L = 5 our Monte Carlo 
value coincides with the exact one (computed by exact numeration of the 2^ configura- 
tions) with an accuracy of less than 10~®. Note also that our MC values converge as the 
size is increased to the exact infinite-lattice value as given by the Onsager solution [Q. At 
T = 8 (high-temperature regime) our representation is less good but still very satisfactory. 
As a function of the size, the gain in computational effort converges and a value of about 
20 is gotten. At the infinite-lattice critical value the results are less spectaculary but still 
of interest. A converged value of about 3 for the gain in efficiency is obtained. At this 
temperature the correlation length for the spin variables diverges and more accurate rep- 
resentations for the solution of Eq.(|l^ are needed. Starting from our basic equation built 
from a transition probability corresponding to local moves we need to resort to approximate 
solutions which contain in some way the collective spin excitations. Alternatively, we can 
change our fundamental equation by resorting to a non-local transition probability density 
and, then, to a new operator H. A natural choice is of course the transition probability of 
the Swendsen-Wang algorithm used here to generate configurations. Preliminary calcula- 
tions show that statistical fiuctuations are indeed strongly decreased. However, to sum up 
analytically all contributions corresponding to the different Swendsen-Wang clusters (action 
of H on ip) is very time-consuming and the advantages of the method can be lost. Some 
approximate scheme is clearly called for; this is let for future development. Finally, a last 
important point is that the gain in computational effort is found to be systematically greater 
(by about 50%) than the corresponding ratio of variances. This result is a direct consequence 
of the fact that the integrated autocorrelation time known to control the amount of correla- 
tion between successive measurements (see, e.g., 0) has been decreased when passing from 
the bare observable to the renormalized one. Note that a similar behavior has also been 
obtained in applications based on improved estimators 0, 0. Without entering into the 
details, it can be shown that this result is directly related to the fact that the fluctuations 
of the renormalized observable are much smaller than in the bare case. 

The second application illustrates the method in the case of a continuous conflguration 
space (calculation of multi-dimensional integrals). We have calculated a mean energy as 
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it appears in the so-called Variational Monte Carlo (VMC) methods 0. Starting from a 
quantum Hamiltonian Hq (to be distinguished from our trial operator H) and a known trial- 
wave function ipT our purpose is to compute the variational energy associated with ipT- 
Ey can be easily rewritten as an average over the probability distribution ipT^^ Ey =< El > 
where El = Hqi/jt/'^Pt is called the local energy. Here, we consider the case of the Helium 

atom described by the Hamiltonian Hq = — l/2(Vi + V2 ) — 2/ri — 2/r2 + l/ri2 (atomic 
units) with usual notations. As trial wave function a standard form has been chosen 

ipT{ri,f2) = exp[ '^^^^ c(ri + r2)]ls(ri)ls(r2) (15) 

i + ori2 

where ls(r) is the Hartree-Fock orbital as given by Clementi and Roetti [§] and the varia- 
tional parameters have been chosen to be a = 0.5,6 = 0.522, and c = 0.0706. As already 
remarked a natural choice for the trial Hamiltonian H is a Schroedinger operator admitting 
'0T as ground-state, Eqs. ([To|JTl|) . Regarding ip we have chosen a form similar to the trial 
wavefunction multiplied by some function of the potential energy. Configurations are gener- 
ated using a standard Metropolis algorithm with local moves constructed using a Langevin 
equation [0 . Results are presented in Table It is seen that the introduction of the renor- 
malized local energy increases the efficiency of the Monte Carlo calculation by about one 
order of magnitude. 

In the last application it is shown that the method can even be used in exact (zero- 
temperature) quantum Monte Carlo (QMC) calculations. In QMC a combination of diffusion 
and branching process is used to construct a stationary density proportional to V'r'^o where 
ipo is the exact unknown ground-state wavefunction. By averaging the local energy over 
this distribution an estimate of the exact energy Eq is obtained 0. Although the analytical 
form of the stationary density is no longer known, a renormalized function whose average 
is identical to that of the bare local energy can still be defined. El = El + {H — Eq)iP / %Ijt 
where H admits ipx as eigenvector, Hipx = 0. Calculations have been done using the exact 
Green's function Monte Carlo of Ceperley and Alder |p . Results are presented in Table 0. 
They are of a quality similar to that obtained in the variational case. About one order of 
magnitude in computer time has been gained. 

To conclude we have presented a simple and powerful method to greatly increase at little 
cost the efficiency of Monte Carlo calculations. The examples presented have been chosen to 
illustrate the great versatility of the method (discrete and continuous configuration spaces, 
classical or quantum Monte Carlo, local or non-local Monte Carlo updates). Although our 
examples have only been concerned with total energies, let us emphasize that the zero- 
variance principle is valid for any type of observable including important quantities such as 
local properties other than energy, differences of energies, spatial correlation functions, etc... 
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TABLES 

TABLE L Internal Energies for the 2D-Ising model at different temperatures. N number of 
sites. Statistical uncertainties on the last digit are indicated in parentheses. 

Si^i 5 x 5 10 X 10 20 x 20 25 x 25 oo x oo 



a{Ey/N 
a{Ef/N 
Ratio of variances 

< E/N > 

< E/N > 

Gain in computational effort" 

< E/N > Exact value 



1.789(1) 
0.0125(4) 
~ 143 



1.777(2) 
0.0061(1) 
~ 291 



1.78(1) 
0.0060(2) 
~ 297 



1.79(1) 
0.0061(2) 
~ 293 



-3.902044(31) -3.902200(55) -3.90217(21.4) -3.90242(29) 
-3.902020(2.4) -3.902229(3) -3.90225(1.2) -3.90222(1.5) 

~ 167 ~ 336 ~ 318 ~ 360 

-3.9020214... 



-3. 9022331. ..'^ 



Tr = 4.53837... 



a{Ey/N 18.581(4) 25.97(3) 33.1(2) 

a{Ef/N 0.215(2) 4.85(1) 16.5(1) 

Ratio of variances ~ 86 ~ 5.4 ~ 2.0 

< E/N > -3.07334(13) -2.95214(33) -2.8902(12) 

< E/N > -3.07345(1.3) -2.95236(13) -2.8908(7) 
Gain in computational effort" ~ 100 ^ 6.8 ~ 3. 

< E/N > Exact value -3.0734396... 
r = 8 



35.3(2) 
16.9(2) 
~ 2.1 
-2.8800(14) 
-2.8788(8) 
~ 3.1 



-2. 8284271... 



a{Ey/N 13.17(1) 10.96 

a{Ef/N 0.041 0.455 

Ratio of variances ~ 321.2 ~ 24. 

< E/N > -1.16440(33) -1.11556(48) 

< E/N > -1.16348(1.6) -1.11502(8.2) 
Gain in computational effort" ~ 425 ~ 34 

< E/N > Exact value -1.1634926... 



11.1(2) 

0.8(1) 

~ 13.9 
-1.1156(20) 
-1.1145(4.4) 

~ 20.7 



10.9(3) 

0.9(1) 
~ 12 
-1.1165(25.6) 
-1.1145(5.6) 

~ 20.9 



-1.1145444...'' 



"■ See text for definition. 
''Ref. [|. 
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TABLE II. Energy of the Helium atom. All quantities are given in atomic units. Statistical 
uncertainties on the last digit are indicated in parentheses. 

Variational Monte Carlo 



Ratio of variances 

<El> 

< El > 

Gain in computational effort" 



Exact Green's Function Monte Carlo 



a{ELf 

Ratio of variances 

<El> 

< El > 

Gain in computational effort" 
Exact energy 



" See text for definition. 
''Reference Pi 

CT 



^Reference [10| 



0.0409(2) 
0.00688 
~ 5.9 
-2.89671(4.8) 
-2.89674(1.6) 
~ 9 



0.0411(9) 
0.00855(8) 
~ 4.9 
-2.903745(99) 
-2.903734(33) 
~ 9 

-2.903724377... 
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